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Do liquid drops roll or slide on inclined surfaces? 
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A solid sphere is likely to roll, while a rectangular box is likely to slide, on an inclined surface. 
In contrast, a liquid drop moving on an inclined surface can exhibit a variety of shapes and hence 
complex but interesting dynamics. Combining lattice Boltzmann method for hydrodynamics with a 
diffuse interface model, a hybrid numerical scheme is used to study the dynamics of binary fluids on 
an inclined plate under the action of gravity. Using a triple decomposition of the velocity gradient 
tensor inside a drop, the vorticity associated with the rolling motion is distinguished from that of 
shearing motion. The average angular velocity based on this residual vorticity becomes significant 
when the external fluid viscosity is reduced and the shape of the drop approaches a circle. For 
a given slip length and viscosity ratio with the outer fluid, irrespective of the equilibrium contact 
angle, plate inclination or Bond number, a universal curve is observed for the amount of rotation 
as a function of the drop shape characterized by the isoperimetric quotient. Results from a large 
number of simulations over a wide range of parameters are shown to lie on this curve. The rolling 
motion is also found to be strongly dependent on the slip length at the contact line. 

PACS numbers: Valid PACS appear here 
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I. INTRODUCTION 

A solid sphere on an inclined surface rolls down under 
the action of gravity. A rectangular object however is 
more likely to slide down. The choice of sliding versus 
rolling motion is determined by the shape and weight of 
the body and the frictional forces at the supporting sur- 
face. On the other hand, a liquid drop of a given volume 
on an inclined surface can attain a variety of shapes. The 
static shape depends on the equilibrium contact angle, 
solid surface characteristics including the pinning loca- 
tion, interfacial tension between the substrates, gravity 
and plate inclination, while a moving drop has to contend 
in addition with viscous and inertial stresses More- 
over, the reaction forces and moments provided by the 
supporting surface are distributed and depend strongly 
on the shape, see e.g. @, These differences between a 
solid and a liquid make the study of the latter complex, 
but pose an interesting question, viz. whether a liquid 
drop sitting on an inclined solid surface will roll, slide, or 
do both. 

The motion of the liquid drops on solid surfaces is 
much studied, especially in the two limiting cases of per- 
fect wetting and zero wetting. When the contact angle is 
small, the lubrication approximation of the Navier-Stokes 
equation is a good model to describe the dynamics, 
with a prescribed parabolic velocity profile. The sliding 
motion of the drops and the associated instabilities of 
the receding front has received much attention theoret- 
ically and experimentally, see e.g. [1, Q- At the other 
limit, when the contact angle approaches 180°, a com- 
plete rolling of the drop is observed under certain situa- 
tions An anomalous increase in the speed of smaller 
drops has been observed on super hydrophobic surfaces 
3 and explained based on the scaling arguments in the 
model of 0. The shape of these droplets as they roll 
down has been studied by While these two limits 



are well studied, intermediate contact angles are studied 
much less. They are harder to analyze since they do not 
lend themselves to simplifying approximations. While 
analytical solutions arc practically impossible, numerical 
solutions pose considerable challenges due to the multiple 
length scales present and the coupling between the evolv- 
ing field and the interface shape Here we analyze 
the entire spectrum of shapes for two-dimensional drops, 
for a range of the relevant non-dimensional parameters. 
A splitting of the motion into sliding, shear and rolling 
leads to a better understanding of the dynamics. 

The competition between rolling and sliding motion of 
droplets on hydrophobic surfaces has been investigated 
experimentally [l2f, ■ In these studies the effect of sur- 
face coating and roughness on the internal fluidity of the 
droplet was of primary concern. The velocity inside the 
droplet was obtained by particle image velocimetry, and 
the contribution to slip versus to roll was evaluated for 
an accelerating drop on an inclined surface. Rolling mo- 
tion is clearly observed in molecular dynamics [T^ . and 
lattice Boltzmann, simulations [ist on cylindrical drops. 
However, such studies have concentrated on the total ve- 
locity of the droplet and its dependence on the driving 
force and the contact angle. With a somewhat different 
goal, we investigate here the motion of cylindrical drops 
of various contact angles on smooth surfaces, where the 
effect of fluid properties like density and viscosity, and 
also of slip are evaluated. 

The idea of splitting the motion into shear and roll 
is standard in fluid mechanics. However, the standard 
splitting docs not distinguish global rotation from local 
rotation of the fluid element. Such distinction is neces- 
sary to understand the global dynamics of a drop, so we 
distinguish between these, in the manner introduced by 
[l6| in a different context. A splitting of the velocity field 
in a drop into slip and roll was done recently by [ol , al- 
though in a different manner than done here. The linear 
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part of the velocity profile on a line passing through the 
center of the drop was attributed to rotation whereas we 
will see in section IIIDI that this would overestimate the 
rotation inside the drop. In addition, the study of [l7| is 
valid for drop sizes smaller than the capillary length, so 
large deformations from the circular shape were not un- 
der consideration, while the present approach has no such 
restriction. Wc show that the shape, and hence the size 
is very important in determining the amount of rolling 
inside the drop. 

Movement of the contact line, by slipping or other- 
wise, is imperative for a moving drop unless the contact 
angle is exactly 180° @. Though the exact mechanism 
of contact line movement is not understood [l^, it is 
generally believed that the macroscopic behavior of the 
drop is independent of the assumptions at the contact 
line. However, it is clear that even a small amount of 
global rotation can manifest itself as tank-trcading [l^ 
near the contact line, which means that contact line be- 
havior is not local in nature, and must be consistent with 
the macroscopic motion. We demonstrate this using the 
diffuse interface (DI) model. Most earlier studies have 
concentrated on the rolling motion near the contact line 
(20I [2l| while we look at the rolling motion in the bulk 
of the drop. The association between the two, if any, 
implies the nonlocal hydrodynamic effects of the contact 
line movement 

Understanding the kinematics is not just a curiosity 
but has practical relevance too. For example, rolling 
droplets play an important role in self-cleaning devices. 
As they roll, they pick up and remove dirt as observed on 
hydrophobic surfaces • Hence it is desirable to know 
under what situations one can maximize the rolling mo- 
tion inside the drop. 



II. THEORY: DIFFUSE INTERFACE MODEL 
AND HYDRODYNAMICS 

We use a coupled system of equations describing the 
hydrodynamics of a conserved order parameter ijj and the 
conserved momentum density pu, where p and u are the 
total density and the local fluid velocity. 



A. Landau-Ginzburg theory 

For a binary fluid system consisting of species / and // 
with local densities nj and n//, the order parameter is 
defined as the normalized density difference, ip = ^^^^ 



nii+ni 

which quantifies the local composition. The equilibrium 
thermodynamics of the fluid is described by the Landau 
free-energy functional [H, [131 



(./W + ||V^P)dr, 



(1) 



and is approximated as = ^'4'^ + ^4''^ with A < Q 
and B > Q. The three parameters A, B, and K con- 
trol the interfacial thickness and interfacial energy of the 
mixture. The second term of Eq. [T] involving the square 
gradient gives a free energy cost to any variation in the 
order parameter, and is related to the interfacial tension 
between the two fluid phases [2^ . Two uniform solutions 
ijj = A/B can coexist across a fluid interface. For a 
planar interface; the concentration profile between the 
two bulk phases is given by 



^(z) 




(2) 



where z is the coordinate normal to the interface while 
^ = ^ ^ determines the interfacial thickness. The 
energy associated with this profile in excess of the en- 
ergy in the bulk, defined per unit area, provides the 

interfacial tension 7 = The corresponding 

chemical potential is given by the variational derivative 
of the free energy with respect to the order parameter 
fi = SF/S-tp = Alp + Bip^ - KV'^ip. Gradients in the or- 
der parameter produce additional stresses, which follow 
from the relation tj)^ fx = V • cr"^ [2^, including Laplace 
and Marangoni stresses due to a fiuid-fluid interface. 



B. Governing equations 

The order parameter is described by a Cahn-Hilliard 
equation (CHE), which includes advection by fluid flow 
and relaxation due to chemical potential gradients. 



dtij + V-iuij) = V-(MV/i). 



(3) 



where r stands for the spatial dimensions. The first term 
represents the local free energy density of the bulk fiuid. 



The mobility M is the constant of proportionality in the 
linear phenomenological law relating the thermodynamic 
flux of tp to the thermodynamic force V/x. The order pa- 
rameter dynamics is coupled to a Navier- Stokes equation 
(NSE) [27| with additional stress densities arising from 
the order parameter. For an incompressible fiuid, the 
dynamics is governed by 

dtipu) + V- (puu) = -Vp + r/V^u + pVp, + G (4) 

together with the continuity equation for the density. In 
the above, p stands for the isotropic contribution of the 
pressure, 77 is the shear viscosity and G is the gravita- 
tional force density. 



C. Numerical Algorithm 

We now briefiy review the numerical algorithm used to 
solve these coupled equations. We use a hybrid algorithm 
by combining the lattice Boltzmann method for hydro- 
dynamics and method of lines for the order parameter 
dynamics. The interested reader is referred to [l^ for a 
detailed description. 
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The Lattice Boltzmann (LB) method for solving the 
Navier-Stokes equations is modified to include force den- 
sities such as the divergences of order parameter stresses 
and gravity [1^. In a standard DdQn LB model where 
the velocity space is discretized into n components in d 
dimensional space, the discrete form of the Boltzmann 
equation reads 



where F(x, i) is an effective force density. The moments 
of the single particle distribution function fi, defined at 
lattice node x with velocity Cj at time t, give the fluid 
mass, momentum and stress densities: 
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where Qiap = CiaCip — clSai}- The collision operator Ljj is 
the discrete form of the collision integral. It controls the 
relaxation of fj to equilibrium, /j*, and may be modeled 
in terms of a single relaxation time t as Lij = Sij/r. The 
viscosity is then obtained as r/ — tc^ where Cg = l/\/3 is 
the sound speed in LB units. 

The spatial discretization of the Cahn-Hilliard equa- 
tion is based on a finite- volume formulation [s^]- For a 
given node, the divergence is written as a sum of fluxes 
defined on the midpoint of the link connecting the node 
to its neighbors. The resulting equation is temporally 
integrated using a Rungc-Kutta algorithm (28l |. 

A liquid drop will typically be of higher viscosity and 
density than its surroundings. Secondly, the contact line 
must be treated with care. In the present work, these 
features are incorporated into the numerical approach as 
discussed below. 



1. Viscosity Ratio 

The differences in properties between the two fluids 
could be introduced at a molecular level, as done by 
[Slllll]- Our approach is macroscopic, and the simplest 
way to introduce a viscosity difference across the fluid 
interface is to prescribe the relaxation time as a func- 
tion of the order parameter. The underlying assumption 
is that the molecular structure of two fluids is the same, 
and is analogous to the introduction of interfacial tension 
using Cahn-Hilliard theory. To test this, we first design 
a model problem, of the laminar pressure-driven flow of 
two fluids in a two-dimensional channel. Fluid I is of 
higher viscosity rji, and occupies the lower portion of the 
channel, while fluid II of lower viscosity rju occupies the 
upper portion. 

In the literature we find two different expressions for 
the relationship between the concentration and the re- 
laxation time. 



FIG. 1: Comparison of 2D channel flow velocity proflles of 
an immiscible binary system and its approximations in the DI 
framework. The thick continuous line is the velocity profile of 
an immiscible binary fluid. In the simulations, the binary fluid 
was approximated as a viscosity stratified fluid flow according 
to Eq. [7] and [H] the results are shown by red stars and blue 
circles respectively. The corresponding analytical solutions 
for viscosity-stratified flow are shown by the red dashed and 
blue solid lines respectively. While LB perfectly reproduces 
the flow for a given viscosity stratification, the use of a thin 
viscosity stratified layer itself is a good approximation for an 
immiscible binary fluid system. 

(i) Effective relaxation time as a polynomial function 
of order parameter, the simplest being linear [H, [s^l. 



r==0.5[T/(l-^)] + [rj7(l + ^)]. 



(7) 



(ii) As Arrhenius suggested, prescribe an effective re- 
laxation time as a product of relaxation times, but raised 
to a power proportional to concentration (ssj 



^11 



(8) 



We test both relationships by simulating the two-fiuid 
fiow described above, first defining the viscosity ratio 
simply by rjr = rji/rjjj — tj/th. The velocity profile 
from simulations using the two expressions above are 
compared to the analytical solution for two immiscible 
fluids separated by a sharp interface in Fig. [T] When 
the viscosity contrast is small, both match well with the 
analytical solution. However, when the ratio is large {r]r 
= 10), Eq. [5] is closer to the immiscible result than Eq. 
[T] This is because the effective mixed layer where the 
viscosity varies between ?// and rjn is smaller by the for- 
mer relationship. This is evidenced by the fact that a 
corresponding analytical solution for the laminar veloc- 
ity profile for the parallel fiow of two miscible fiuids with 
a thin mixed region between them agrees in each case 
with the computed result. We have thus shown that this 
approach is a good one for incorporating viscosity con- 
trasts in LB simulations. Note that no ad-hoc fixes are 
needed. We have used Eq. [5] in our calculations. 



2. Gravity 

The hybrid algorithm written using a single particle 
distribution function is modified to include gravitational 
effects as follows. Following [s^, we may write the grav- 
itational force acting on the fluids as, 



G 



(9) 



where Ig is the direction in which gravity is acting. 
Therefore, when = 1, G = pgulg and when ip = —1, 
G = p.g/lg. At the interface when ip = Q, G = 
2(91 ~^ 3//)lg- Thus gi/gii determines the density ratio 
between two fluids. In order to ensure that we are in the 
incompressible limit, we must have \Gz\ << pc^ where 
z is the vertical extent of the simulation domain, which 
means that thermodynamic pressure is large compared 
to the hydrostatic pressure difference. 

In our simulations, wall boundary conditions are ap- 
plied on two sides of the domain. Thus a computa- 
tion of a single component fluid, with gravity prescribed 
along the flow direction, develops a channel flow be- 
tween the walls. For droplet simulations, we implement 
the body force only on one fluid, which is equivalent 
to solving the NSE with the Boussinesq approximation 
as described below. For small density variations, i.e., 
Ap = gj — gii << gj, the Boussinesq approximation 
provides that we may neglect the density variation ev- 
erywhere in the NSE except in the buoyancy term, G of 
Eq. 21 For simplicity we may further absorb the body 
force pgii into the pressure term by redefining pressure. 
Density in the interface region is prescribed as a linear 
function of the order parameter. Note that the validity 
of our simulations is thus limited to situations where the 
Boussinesq approximation holds. 



3. Wetting Boundary Conditions 

The algorithm implemented toget the correct contact 
angle on the wall is based on [33, HI]. The solid- fluid 
surface tensions are introduced by defining the Landau 
free energy functional 



(10) 



where '0s is the value of order parameter at the wall. 
Minimization of this energy functional near the wall gives 
a relation between energy gradient and the gradient of 
the order parameter 



d->ps 



kWtp ■ n. 



(11) 
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where n is normal to the wall. The form 
H^ps is known to be sufficient to produce various wetting 
behavior. By tuning the parameters C and H we can 
modify the properties of the surface. If i7 = we have 
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FIG. 2: Verification of the implementation of wetting con- 
ditions, Eq. 1121 on the walls by comparing the numerically 
obtained contact angle with Eq. 1131 The agreement is good 
except at very large and small contact angles. 



neutral wetting. Nonzero values of H therefore allows an 
asymmetry in the surface value of the order parameter 
and a contact angle different from 90°. Therefore we 
have, 



dfs „ , 
dyjs 



H = KWtl^ ■ n. 



(12) 



It is found sufficient to retain only the linear term of 
the surface energy functional [s^, HI], i.e., to set C = 
0. We use a second order central difference formula to 
calculate the normal derivative of order parameter at the 
wall. Thus, H = K^^^^^, where the subscripts and 1 
represents the 0*'* and first node respectively, which may 
be used to obtain the order parameter at the boundary 
node. The wall is placed at the ^ location, as is usual 
in the bounce back schemes used to represent wall in LB 

procedures [s^. Defining a parameter h = the 
contact angles may be calculated as 



cos 9 



(1 + /l)3/2 _ (1 _ /i) 



3/2 



(13) 



In addition po ~ pi is imposed to ensure no order pa- 
rameter flux into the wall, up to second order accuracy. 
Also, advection terms have been carefully discretized to 
preserve mass conservation upto machine accuracy psj . 
The implementation of the wetting properties of the wall 
is verified by comparing the static equilibrium angle ob- 
tained from the simulation to that from Eq. [13] as shown 
in Fig. [2] Deviations are seen at very large and very small 
contact angles and this appears to be consistent with the 
literature [33, HI] (both used kinetic schemes for CHE) . 
Both the grid size and the interfacial thickness were dou- 
bled separately to ensure that the deviations seen are not 
due to the contact line pinning. One reason for the de- 
viation may just be the error involved in the derivation 
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FIG. 3: (a) tp, order parameter field in tiie entire simulation 
box, two colors for the two different fluids with a thick inter- 
face, (b) Vcm, velocity field in the center of mass frame of the 
drop, the rolling motion of the drop may be observed, (c) uj, 
the vorticity field inside the drop, regions concentrated vor- 
ticity may be noticed near the contact line, (d) W, a measure 
for the Weiss criterion to distinguish regions of high shear 
from vortical regions, (e) m, the kinematic vorticity number, 
again a measure of highly vortical regions, (f) cures, residual 
vorticity showing the regions associated with solid body rota- 
tion. 



of Eq. [13] where h is assumed to be small and Taylor se- 
ries expansions are used. In spite of this, we restrict our 
simulations as far as possible to the wide range of inter- 
mediate contact angles where we are sure the simulations 
capture the relevant physics in this respect. 



D. Measure of Rolling Motion 

Before discussing the simulations and results, we take 
a typical simulated drop and discuss how the slide, shear 
and roll may be estimated. The droplet is illustrated in 
The rolling motion inside the drop is evident in 



3(a) 



the corresponding velocity field, plotted in the center of 
mass reference frame of the moving drop as illustrated in 
|3(b)[ Our objective is to quantify this rotation. The first 
quantity to look at would be the vorticity field, shown in 
|3(c)[ The velocity gradient tensor is usually split into a 
symmetric part S and an antisymmetric part CI as 



Vu = S + r2. 



(14) 



S(a;, y) for shear tensor is a measure of the deformation of 
a small fluid element located at {x,y), while the vorticity 
tensor fl{x, y) identifies the angular velocity of the fluid 
element, with {x,y) as the center of rotation. There are 
different measures available in the literature to identify 
regions of high vorticity relative to shear (ioj , and most 
of them have been derived in the context of turbulence. 
A commonly used measure is W = — ||S| 

Weiss criterion. This is illustrated in Fig |3(d) 



2 of the 
for the 



drop under consideration. Another measure is the kine- 
matic vorticity number, defined as m = ||r2||/||S|| where 
= tracc([.] • [.]^)^/^. This quantity is plotted in Fig 
Different criteria have been developed later on to 



3(e) 



dimensions to the Weiss criterion [411 . However the main 
drawback of these measures is that they all estimate vor- 
ticity, which does not in general give a direct measure of 
rolling motion. This is because the vorticity a; = V x u is 
a local quantity which includes both solid body rotation 
and shearing motion of a fluid element, and thus cannot 
distinguish between them. The residual vorticity which 
we describe below is shown in [3(f)[ Although all three 
measures broadly describe the region of high rotational- 
ity in a similar fashion, only the last is good for obtaining 
a quantitative estimate of solid body rotation. 

To demonstrate that these conventional measures of 
vorticity are inadequate in describing the global rolling 
motion inside a drop, let us discuss a simple flow con- 
figuration, namely, a shear flow with u = jy where x is 
the flow direction and y is the gradient direction. The 
velocity gradient tensor Vu splits into its symmetric and 
antisymmetric parts as follows: 



Vu = 



7 




7/2 
7/2 





7/2 



-7/2 




= s + o. 



Although n is non-zero, there is no rolling motion, which 
tells us that we need a modified vorticity to characterize 
rolling motion. 

To demarcate a coherent vortex correctly from high 
shear regions, Kolar [l6j proposed a scheme to remove 
the 'shear' vorticity from the total vorticity. He proposed 
a triple decomposition of the relative motion of a fluid 
element, where the velocity gradient tensor is split into a 
straining part, a rigid body rotation and a simple shear 
flow part. For clarity, we illustrate these pictorially in 
Fig. [?] Given a velocity field in two dimensions, u, the 
velocity gradient tensor (Fig. 4(a) ) is a 2 x 2 matrix 







-{ 


\Vx 


^y) 





Ux 







.(15) 



The symmetric part can be diagonalized to give 
("^o^ ^5/2) ■^here s = ^4u^ -|- {uy + VxY- This is the 
strain rate tensor in the principal coordinates, see Fig. 
[4(b)[ which represents the total straining of the fluid ele- 
ment. The rotation tensor being antisymmetric will not 
change with a rotation of the coordinate system, and re- 
mains as (^^^2 0^^) ^'^'^'"c 1^ — Ux ~ Uy^ the vorticity. 
Therefore, in the principal axis coordinates, the velocity 
gradient tensor is, (^^/2 -s/2 ) ■ '^'^^ rotate the coor- 



dinate system further by 7r/4 (Fig. 4(b)). This frame is 
called a basic frame of reference (BFR), and the velocity 
gradient tensor in this frame is 

(s-w)/2\ /O s/2\ / -cj/2 
{s + uo)/2 j ~ i^s/2 )^\uj/2 

In this reference frame, the contribution due to shear is 
maximized in a triple decomposition of Vu. One may 
write [llj 

Vu = Vu, 



define a vortex exactly and many of them reduce in two 
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Original fiow, BFR 





(c) Example I 



Originai fiow, BFR 





(d) Example II 

FIG. 4: (a) Standard velocity gradient decomposition into 
symmetric and antisymmetric parts for a a flow with Ux =0.1, 
Uy — 0.3, Vx ~ —0.2, Vy — —0.1. (b) The strain rate tensor in 
the principal coordinate system is rotated by 45° and added to 
the antisymmetric tensor to generate the same flow in BFR. 

(c) The same flow is decomposed into a simple shear flow 
and residual flow, the residual being purely rotational flow. 

(d) A strain dominated flow, = 0.1, Uy — 0.2, — 0.1, 
Vy = —0.1 is decomposed into simple shear and residual flow, 
with residual consists of only straining flow. 



Only one of the residual terms Sresiduai or flresiduai will 
be non-zero. The residual matrices above are constructed 
as explained below. The residual vorticity and strain may 
be written respectively as 

UJres = if \s\ > \uj\ 

- sgnH [|c.| - |5|] if \s\ < \ul 

Sres = Sgn(s) [\s\ - \uj\] if |s| > |w| 

= if |.s| < \uj\. 
An example flow which is vorticity dominated, where 



|s| < |a;| is illustrated in figure 4(c) For this case we 
may write 





s/2 , 








sgn{uj)\§\ 


2 



sgn(ui)— — — 



-sgn{uj)\i 




/ \ 111)! — Isl 
-sgn[ijj) 2 





where the first two matrices will combine to produce a 
simple shear flow. The residual straining, shown by the 
third matrix, is zero. The portion of the velocity gradient 
tensor contributing to solid body rotation is given by the 
fourth matrix as illustrated in Fig. 



4(c) 



On the other 

hand, for a flow which is strain dominated, i.e. \s\ > \uj\, 
an example of which is demonstrated in Fig. |4(d)[ we 
write 





sgn{s)\^. 







sgn{s)\ 


sgn{s) 




sgn{s) 



Again the first two matrices on the right hand side to- 
gether produce a simple shear flow. The remainder is 
seen in the third matrix to be purely straining, as illus- 
trated in Fig. |4(d)[ We may now see that for the simple 
shear flow discussed above, \s\ = \uj\ and uJres = 0, while 
for a solid body rotation, lo = ujres- If |s| < the 
flow is vorticity dominated, and the residual tensor will 
consist of only rotation, and vice versa. Therefore resid- 
ual vorticity can characterize the rolling motion inside a 
drop. 

This separation of the shear vorticity from residual is 
different from the shear and curvature vorticity used by 
the atmospheric science community. There the shear vor- 
ticity is defined as — where u is the local velocity along 
the streamline and n is normal to the stream line. The 
remainder is defined as the 'curvature vorticity', which 
is associated with the curvature of streamlines. Note 
that these values will not be Galilean invariant. This 
procedure, when applied to a solid body rotation, pre- 
dicts equal values for both shear and curvature vorticity, 
though there is no shear component present in solid body 
rotation. Hence, we will not benefit here from this pro- 
cedure. The residual vorticity that we use does not suffer 
from this drawback. 

In contrast, in an irrotational vortex where Ur = and 
ug = l/r the residual vorticity is zero because vorticity 
associated with shear and rotation are equal and of op- 
posite signs. However the curvature vorticity is non-zero 
showing the swirling motion of ffuid elements. We may 
neglect such a contribution as we do not expect a point 
vortex like motion inside the drop. In other words, in the 
strict limit of Stokes flow, the velocity field inside a drop 
is constituted by growing harmonics alone. Therefore it 
is reasonable to consider that the residual vorticity gives 
the correct quantitative measure of rolling motion inside 
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a drop. Note that neither residual vorticity nor curvature 
vorticity are complete measures of rolling motion. One 
may use residual vorticity in low Re flows and curvature 
vorticity in high Re flows. 

By splitting the vorticity into two parts, one can iden- 
tify the regions of shear. This could have also been iden- 
tified by looking at the shear rate or viscous dissipation. 
However, such an approach will miss a very important 
factor to the motion which is solid body rotation which 
will not produce shear, but is important in determining 
the dynamics. Hence looking at the actual and residual 
vorticity plots gives an idea about overall motion, strain- 
ing regions and rotating regions. As we will see that 
solid body rotation is indeed an important ingredient in 
the dynamics. 



III. RESULTS AND DISCUSSION 

Simulations have been performed in a box of dimen- 
sions 512 X 256 X 1 LB units for a cylindrical drop. Wall 
boundary conditions are applied on two sides and pe- 
riodic boundary conditions are applied on the other two 
sides. The simulation is initiated with a semicircular drop 
of radius L = 60 sitting on one wall, which is inclined at 
an angle a to the horizontal. In response to gravity and 
surface forces, the drop starts moving on the solid sur- 
face. We impose a smooth surface which thus does not 
display any hysteresis. The simulation is continued till 
the drop reaches a steady state velocity V. Before com- 
piling all the results into a unified framework, we first 
examine the effect of varying one physical quantity at a 
time. We define the Bond number as Bo = L'^\G\/(j, 
the Reynolds number as Re = LVp/rj and the Capillary 
number as Ca = rfV/a. 

The effect of increasing gravity on the steady state 
drop shape and streamline patterns, total and residual 
vorticity and angular velocity are illustrated in Fig. [S] 
A larger driving force means larger deformations, so the 
drop deviates further from its equilibrium shape at zero 
plate inclination. The streamlines are plotted in the cen- 
ter of mass frame. Fixing the center as the innermost 
streamline, vorticity and residual vorticity along different 
streamlines are plotted as functions of azimuthal angle, 
thus mapping the entire vorticity field inside the drop. 
The thick interface and the contact line region are ex- 
cluded from quantitative consideration as there may be 
spurious velocities generated due to the LB-DI model 
[4^ . In Fig. 5(a) the total vorticity is seen to lie within 



a small range almost everywhere in the drop, but as dis- 
cussed above, we may not use this to determine whether 
there is a solid body rotation. In this case, the resid- 
ual vorticity indicates that the bulk of the drop is indeed 
in solid body rotation. The total and residual vorticity 
values are comparable, showing that there is hardly any 
shear vorticity. However, in a region near the rear con- 
tact line, shear vorticity dominates. Thus the up-down 
symmetry is broken. Finally the angular velocity based 



on residual vorticity, Vres = fiOres where r is taken as the 
radial distance from the center of the innermost stream- 
line, is plotted as a function of azimuthal angle. A perfect 
solid body rotation, such as that of a solid wheel would 
have appeared as concentric circles in this plot. At small 
Bond number, we do have something like this, except for 
a slight loss of symmetry in that the outer streamlines 
move faster at the top and slower at the bottom. In the 
case of large Bo, the drop is elongated normal to grav- 
ity, with a clear breakdown in left-right symmetry in its 
shape, as illustrated in Fig. 5(b) Except for the very 



center of the drop there is no resemblance to solid body 
rotation or even to tank-treading. This is reflected in the 
the angular velocity plot as well. The residual vorticity 
is higher in the direction of elongation. Interestingly the 
residual vorticity is now higher near the rear of the drop, 
exactly where it was lower at low Bo. This is because 
at higher gravity the rear of the drop has a tendency to 
lift off the surface. A given fluid element accelerates and 
decelerates significantly as it moves on a streamline. 

If a solid body is rolling on an inclined surface with an 
angular velocity of A^. then the corresponding vorticity is 
2N . Therefore, we can find the average residual vortic- 
ity inside a drop and calculate a corresponding forward 
velocity of the drop corresponding to the roll as 



Vroll 



Average(wres) Height of the drop 



ing 



(16) 



Here we take the radius of the drop to be half of the 
maximum height. Then a quantity called percentage ro- 
tation, denoted by %R. is calculated based on the total 
translational velocity V of the drop, as 



%R = VrolHng/V X 100. 



(17) 



This is different from the roll versus slip velocity as de- 
fined in [l3| , where the velocity profile at a single location 
is considered and the definition does not distinguish the 
shear vorticity from residual vorticity. In Fig. [5]one may 
see that increasing gravity increases the translational ve- 
locity by an order of magnitude, as seen in the increase 
in Reynolds number, but the associated deformation re- 
duces the percentage rotation %R by 10%. 

Therefore one may infer that an important property 
that determines the motion of a drop is its geometrical 
characteristics. Needless to say, the geometry is in turn 
determined by the volume, the contact angle, the grav- 
itational force and the plate inclination, apart from the 
viscosity and density ratios. The present study is thus 
valid over a wide range of parameters in contrast to Q 
wherein a spherical drop deformed by incrcnicntal grav- 
ity was studied. The crucial assumption in |9| was that 
the deviation of the shape from a sphere is very small. 
Relevant length scales of the deformation as a response 
to gravity were thence derived. These scaling arguments 
break down when 9e 7^ 180° due to a finite contact area 
as we have in our simulations. Also, since we do not re- 
strict our analysis to small Bo, the changes in the surface 
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FIG. 5: Effect of gravity on the drop shape, streamhne patterns, vorticity (oj), residual vorticity (ujrea) and residual angular 
velocity (vres) are illustrated in a coordinate frame moving with the center of mass of the drop. A streamline of a given color 
in the left-most figure is shown in the same color in the three polar plots of lj, Ldres and Vres- The azimuthal angle, measured 
from a line parallel to the solid plate, corresponds to that of the streamline, and the radial location at a given polar angle 
indicates the magnitude of the respective quantities. The magenta lines represent -0 — 0.9, ifj — and '0 = —0.9, showing the 
thick interface. The drop is moving on a surface inclined to the horizontal, and the black dashed line indicates the direction of 
gravity. The red dashed line is normal to it. 



energy need not scale with that of gravitational potential 
energy unlike in 

We now analyze the effect of contact angle alone as 
illustrated in Fig. [S] for drops of same volume. Here 
gravity is adjusted so that the drop attains the same ter- 
minal settling velocity in all cases and hence the same 
Reynolds number. We thus ensure that effects of in- 
ertia are nullified in this comparison. In Fig. [51 one 
may see that no contribution of rotation is present when 
6'e = 42°. Here the entire vorticity of the fluid elements 
can be attributed to that associated with shear. As the 
equilibrium contact angle increases, the percentage rota- 
tion increases with the maximum in the case of an almost 
circular drop. The effect of equilibrium contact angle is 
thus an intuitive result. Compared to the case of Fig. 
|6(b)[ case 6(a) shows a marginal reduction in the per- 
centage rotation despite an increase in the contact angle. 
This can be attributed to the deformation of the drop 
near the wall in the latter case. The total vorticity is a 
function of shape as seen in different cases with a large 



contribution coming from the shear vorticity. 

It is of interest to study how the rolling behavior 
changes as a function of the tilt angle of the plate. This 
is because the ratio of the components of gravity nor- 
mal and tangential to the plate changes. The application 
of the normal component alone does not produce any 
movement of the drop, but both components contribute 
to deciding the shape, and hence the dynamics. The ef- 
fect of plate inclination on the shape, streamlines and 
vorticities and their angular dependencies are illustrated 
in Fig. [T] This illustration is for an equilibrium contact 
angle of 90°. As the plate inclination increases the height 
of the drop increases, and it tends to lift off from the plate 
at some inclination. In turn the percentage rotation in- 
creases, and is highest for a tilt angle of 176°. This is 
another indication that the drop shape is a very impor- 
tant parameter in determining the kinematics inside the 
drop. The presence of corners and deformed parts of the 
drop always increase the shear vorticity locally. 

As the plate inclination changes, not only the ratio of 
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(a) Bo = 0.027, 6»e > 180° (h = -0.7), %R = 17.8, fle : 
0.23, Ca = 0.002 
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(b) Bo = 0.028, 6»e = 138° , %i? = 18.5, Re = 0.23, Ca = 0.002 
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(c) Bo = 0.038, Be = 90° , %R = 8.7, Re = 0.23, Ca = 0.002 

CO 




270 270 

(d) Bo = 0.068, 6»e = 42°,%/? = 1.9, Re = 0.23, Ca = 0.002 

FIG. 6: Effect of contact angle and thence the geometry on 
the roUing behavior. The VoR is larger when the drop shape 
is closer to a circle. In all cases Re and Ca axe kept constant 
by adjusting the Bo and T^r is kept as 10. A larger reduction 
in uires as compared to ui may be observed as 6^ decreases. 



normal to tangential forces ehangcs, but also their mag- 
nitudes. In order to study the effect of this ratio alone, 
the normal force component was artificially varied keep- 
ing the tangential force the same. This corresponds to a 
simultaneous variation in plate inclination and gravity to 
achieve the same settling velocity. The results arc illus- 
trated in Fig. |51 One can clearly see that as the normal 
component of gravity is slowly reduced, the shape be- 
comes more and more elongated in the direction normal 
to the plate and this increases the amount of rotation con- 
siderably. This is yet another indication that the shape 
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(a) a = A,Bo = 0.1, %R = 5.73, Re = 0.57, Ca = 0.005 
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(b) a = 94, Bo = 1.5, %R = 7.3, Re = 8.1, Ca = 0.065 
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(c) a = 176, Bo = 0.1, %R = 10.7, Re = 0.67, Ca = 0.005 

FIG. 7: Effect of plate inclination on the shape and rotation 
behavior of drops. Equilibrium contact angle is 90° and r]r = 
10. A pendant drop is elongated to almost same size as the 
radius, producing more solid body rotation in the drop. 



of the drop plays a big role in the rolling behavior of the 
drop. 

In all these cases, we see that the deviation from a cir- 
cular shape plays an important role in determining the 
dynamics. One can then suitably define a shape param- 
eter to describe the closeness of the shape to a circle, for 
example, the isoperimetric quotient, g, 



47r X Area 
Perimeter^ 



(18) 



This ratio is unity for a circle and is less than this value 
for any other shape, since a circle has the least circum- 
ference for a given area. The percentage rotation in all 
the cases we have computed shows a direct dependence 
on the isoperimetric quotient, as shown in Fig. [9l As 
expected, the percentage rotation is higher for a shape 
which is closer to a circle. It is however of interest to 
note that irrespective of the parameters such as equilib- 
rium contact angle, gravity and plate inclination which 
determines the shape and the deformation, the percent- 
age of roll is a function only of this quantity. The collapse 
of data that we obtain here is a strong indication of this. 
Note that the data in this figure represents results from 
over 100 simulations, spanning a wide range of ^e, ol and 
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(b) a = 4, So = 0.53, %R = 14.5, Re = 4.2, Ca = 0.03, Normal 
component of gravity is I/IO**" of the above. 

FIG. 8: Effect of normal component of gravity on the shape 
and rolling behavior of drops. Equilibrium angle is 90° and 
rir = 10. The tilt angle is chosen as 4°. To differentiate the 
effect of the tangential component of gravity, the normal com- 
ponent of gravity in (b) is artificially suppressed to 1/10*'* of 
its value. However, the tangential component being main- 
tained the same in both (a) and (b) yields a similar settling 
velocity. The percentage rotation can clearly be very different 
even at the same settling velocity, due to the change in shape. 
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FIG. 9: The variation of percentage rotation with the isoperi- 
metric quotient is illustrated for different sets of simula- 
tions. Each color represents a fixed equilibrium contact an- 
gle. Within each set, plate inclinations vary from a = 4° 
to 176°. Also, Bo ranges from 5 x 10~^ to 1.5 by varying 
gravity and surface tension. The slip length and the viscosity 
ratio are kept fixed. An exponential curve fitted through all 
data points is also shown. Note that moving drops for a wide 
variation in physical properties fall on this curve. 




FIG. 10: Isoperimetric quotient of static shapes compared 
with dynamic drops for the same Bond numbers. The equi- 
librium contact angle is 90° and three different plate inclina- 
tions, 30°, 90° and 135°, are chosen for comparison. Here 'FP' 
stands for 'front pinned', 'RP' stands for 'rear pinned' and 'D' 
stands for dynamic cases. The shape factor of a dynamic drop 
lies in between those corresponding to front pinned and back 
pinned static shapes 01 . This behavior breaks down at large 
Bond numbers where inertia is higher. In that case, the mov- 
ing drop is closer to circular than either static shape. 



Bo. However the viscosity, the mobility and the viscosity 
ratio are kept fixed in the simulations shown so far. As 
we will see below, percentage rotation curve gets shifted 
in response to changes in these parameters. We have de- 
fined the outline of the drop as a contour of ?/» = —0.9 
which can be thought of as the inner limit of the inter- 
face. Since we use a combination of LB and DI models, 
there can be spurious interface velocities [i^ and hence 
the data outside this line is not considered. We have also 
tried to calculate this shape parameter from -0 = which 
is theoretically the interface. However this shape fails to 
capture the deformations correctly for large contact an- 
gles. 

It would be interesting to compare the isoperimetric 
quotient of static drops with that of dynamic cases and 
see whether any predictions can be made. This is illus- 
trated in Fig. [TOl As explained in Q we obtain mini- 
mum energy static shapes of drops with either the front 
end or the rear end pinned. As illustrated in Fig. [TUl the 
isoperimetric quotient of the dynamic drops resides be- 
tween that of front pinned and back pinned cases. This is 
interesting because one can make predictions about kine- 
matics inside the drop by the analysis of static drops, 
but only for small Bond numbers. As Bo increases, such 
monotonic variations in the shape factor is violated, ne- 
cessitating the full calculations. 

One of the main drawbacks of the DI models is that 
it imposes a finite thickness of the interface while for 
most macroscopic drops, the interface thickness would 
be negligible compared to any other length scale in the 
problem. In order to ensure that our results are indepen- 
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FIG. 11: A large interfacial thickness is an artifact of the 
DI method. However, here the Ca and %R are shown to be 
independent of Cn. In the above plots the symbols represent 
different sets of simulations; • for Bo — 0.19,6'e — 138°, rjr = 
1, ■ for Bo = 0.04, Oe = 138, = 1, ♦ for Bo = 1.90, 
9e = 90, ?7r = 1 and A for Bo = 0.19, Oe = 138, ?7r = 10. 




0.01 



FIG. 12: Ca is plotted as a function of the nondimensional 
slip length, S. In order to obtain a range of S, the viscosity 
was independently varied by three orders of magnitude and 
mobility by one order of magnitude. The equilibrium contact 
angle is 152°. Larger slip length at the contact line results in 
larger translational velocity of the drop. 



dent of the assumed interfacial width, simulations with a 
range of interface thicknesses were conducted. As shown 
in Fig. [Til both the Capillary number and the percentage 
rotation remain insensitive to interfacial width, to within 
numerical errors. Here Cn = ^/L, the ratio of interfacial 
thickness to the macroscopic length, is the Cahn number. 

Apart from the solid body rotation, which gives a for- 
ward velocity to the entire drop, the contact line moves 
due to the slip provided by the diffusion of the order pa- 
rameter [4^. Balancing the advection and diffusion of 
order parameter across the interface provides a length 
scale for this diffusion as A = ^/ijM. We define a nondi- 
mensional slip length as 5* = X/L. This slip length is 
the same as that used in the slip-induced movement of 
contact line in sharp interface models [i^l, and is not an 
artificial parameter. They actually show that this slip 
length should not be dependent on the interfacial width. 
Hence we can use A as a measure of slip at the contact 




,Bo = 0.38,Ti =1 
I Bo = 0.04, T| 1=10 
► 30 = 0.38,11 1=10 



0.002 



0.004 0.006 

s 



0.008 



0.01 



FIG. 13: Percentage roation, %R, is plotted as a function of 
non-dimensionalised slip length, S. The simulation parameters 
are same as in Fig. 1121 A strong dependence of the rolling 
behavior on the slip length may be observed. 



line. This means that either mobility or viscosity can 
be independently or simultaneously varied to change the 
slip at the contact line. Slip length is here defined using 
the viscosity of the drop rjj. In principle, the external 
fluid plays a very important role and an effective viscos- 
ity to define the slip length may need to account for the 
viscosity ratio rjr- For example a geometric mean of the 
two is used in |44{ . However, we refrain from using this 
relation as it lacks a physical significance. 

Both viscosity and mobility are independently varied 
by at least one order of magnitude in Fig. [T^to obtain 
a range of S*. As the slip length increases, the Ca also 
increases as illustrated in Fig. 1121 Intuitively, a slipping 
drop on an inclined surface will roll less. This is veri- 
fied in our simulations as shown in Fig. [13] wherein the 
importance of slip length in determining the amount of 
rotation inside the drop may be inferred. And this de- 
pendence appears to be exponential. Larger percentage 
rotations than those shown, which would correspond to 
smaller slip lengths could not be obtained reliably with 
the present numerical simulations. 

In our simulations, the nondimensional slip length, S, 
varies from 10^'^ to 10^^. In the light of experimental 
evidence where slip length varies from nm to /xm (45j 
we expect that our observations remain valid for a range 
of drop sizes. For macroscopic drops smaller than the 
capillary length, this ratio is very small and hence a 
larger fraction of rolling motion may be expected than 
those seen here. It is worth mentioning however that 
slip lengths of 10 — 100 of micron have been reported on 
patterned surfaces [i^ or when lubricating gas layers are 
present [i^l or on super- hydrophobic surfaces 22 1. Since 
we concentrate on the bulk motion of fluid elements, we 
expect that our simulations are relevant in several prac- 
tical applications. 

It is interesting to note that viscosity does not have 
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FIG. 14: Though viscosity affects the overall dynamics, the 
contribution of rolling motion to the dynamics remains the 
same. Here a — 30°, 6'^ = 152° and r]r ~ 10. Since viscosity 
plays a direct roll in the slip of the contact line, mobility was 
simultaneously adjusted to obtain same slip length to isolate 
the effect of viscosity alone. 



an explicit dependence on the percentage contribution of 
rolling in the total motion. Viscosity enters the prob- 
lem in two ways, via the Reynolds number, and via the 
slip at the contact line. The former effect is illustrated 
in Fig. 1141 In order to effect changes in viscosity in 
the stress term alone, the slip length was maintained the 
same in the two cases shown by simultaneously changing 
the mobility by the appropriate factor. As a consequence 
of change in dissipation, the settling velocity is different. 
But the percentage contribution of rolling is the same in 
the two cases showing that the viscosity does not explic- 
itly affect the rolling motion. 

Neglecting inertial terms in NSE, one may balance the 
gravitational driving forces to viscous dissipation to ob- 
tain [i,|4|, 



(19) 



where Ag = cos 0^ — c:os 0a , which is small for most of our 
simulations. A plot of Ca vs Bo is shown in Fig. [TSlwhich 
corresponds to the data described in Fig. [S] This scaling 
holds good for all our simulations. The small deviations 
seen are due to inertial effects not accounted for in Eq. 
I19[ as well as by the intercept, as shown below. Having 
seen the effect of contact-line slip on the drop dynamics, 
we now make a connection between the slip velocity and 
the above scaling relation. For this we use the generalized 
Navier's boundary condition [Tsl. l49j. 



—— = [dnu] H d^ip 



(20) 
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FIG. 15: The linear relationship between Ca and Bo is il- 
lustrated for the data given in Fig [9] This data consists of 
different contact angles, plate inclinations, gravity and surface 
tension values. 
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FIG. 16: Both Capillary number, Ca and modified Capillary 
number, Cau are plotted against Bo. The relative standard 
deviation (RSD) is small for CaM justifying Eq. 1221 Both 
viscosity and mobility are varied by at least one order of mag- 
nitude to produce range of slip lengths in these simulations. 
Equilibrium contact angle is 152° and Cn = 0.07. Value of 
/3 is adjusted in Eq. [22] to obtain the least RSD, hence veri- 
fying the role of slip length. Since, 0.1 < Re < 70 for these 
simulations, a wider distribution at large Bo may be related 
to the unaccounted inertial effects. 



to make quantitative estimates of the sliding velocity. In 
the above equation n is a direction normal to the wall 
and C{i(j) ~ Kdni-' + da^f/dip, K is the coefficient in 
the free energy functional and a^f is the intcrfacial free 
energy per unit area at the fluid-solid boundary. The 
second term is called the uncompensated Young's stress 
and one may see that 



dx[C{tp)dx'4>\ = ^{cosOd — cos^e), (21) 



where int means 'across the interface' and 9d is the dy- 
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FIG. 17: Effect of surrounding fluid viscosity on Ca. As the 
viscosity ratio r]r increases, which corresponds to a reduction 
in the viscosity of the surrounding fluid, the drop translates 
faster. 



namic contact angle. We can calculate an order of mag- 
nitude estimate of from the above equation as Ca-|. 
Since we are at steady state, both the advancing and re- 
ceding contact lines move with the same velocity as the 
center of mass of the drop. We assume that any varia- 
tion in order parameter and hence the slip is felt over a 
region of interfacial thickness ^ while the associated slip 
length A is the same slip length calculated in a DI model. 
Accounting for a possible pre-factor in the addition of 
scaling estimates in Eq. [19] we obtain, 



Ca 



Bo, 



(22) 



which explicitly includes the role of slip at the contact 
line. This relationship is verified in Fig. I16[ where the 
change in only the slip length affects the linear relation- 
ship between Ca and Bo. We define Com = [1 + 
By choosing a suitable pre-factor /3 one may see that the 
second term in the above equation explains in part the 
distribution in Ca at a given Bo. In other words, the rel- 
ative standard deviation, which is the ratio of standard 
deviation to the mean value expressed as a percentage 
comes down dramatically when Eq. [221 is used. 

Finally we investigate the role of the viscosity ratio be- 
tween the drop and the surrounding fluid, in Fig. [TTland 
nsl When the viscosity of the external fluid is reduced, 
the settling velocity and hence the Ca, as expected, in- 
crease. Also the percentage of rolling motion is larger. 
In line with this, one may expect significant amount of 
rolling in case of a water-air system where viscosity con- 
trast is large. As the viscosity of the external fiuid in- 
creases and goes beyond that of the drop, the entire dy- 
namics shifts to the external fluid. The drops slides in 
that case. This too is consistent with intuition, since a 
'bubble' will simply slide in a liquid rather than roll when 
moving on a surface. The changes in the vorticity and 
residual vorticity flelds are illustrated in Fig. [THj One 
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FIG. 18: Effect of external fluid viscosity on rolling behavior 
of the drop. The percentage rotation exhibits a strong depen- 
dence on the viscosity ratio, T^r. As the external fluid viscosity 
comes down, rolling motion inside the drop increases. Low 
viscosity drops in a higher viscosity fluid are seen to almost 
slide on the wall rather than execute a rolling motion. 
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FIG. 19: Change in the velocity and vorticity fields when the 
external fluid viscosity is changed keeping the drop viscosity 
same. Here a = 30° and Be ~ 138°. As the external viscosity 
increases, it start affecting the dynamics more as seen in Fig. 

m 



may observe that, despite the geometry remaining simi- 
lar, the percentage rotation increases when the viscosity 
ratio increases. Therefore the universal curve obtained 
in Fig. [Hi will be shifted appropriately by a change in the 
slip length and viscosity ratio. 
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IV. CONCLUSIONS 

A hybrid simulation method implementing lattice 
Boltzmann algorithm with diffuse interface model is used 
to analyze the drop motion on inclined surfaces under 
gravity. Modification of the model to provide a viscos- 
ity contrast between the fluids, wetting boundary con- 
ditions and to introduce gravity are discussed. By re- 
moving shear vorticity from total vorticity, it is shown 
that residual vorticity is a good measure to characterize 
the rolling motion inside the drop. It is shown that the 
drop shape can be described by a geometrical quantity, 
isopcrimetric quotient, that is primarily responsible for 
determining the fraction of forward motion accruing from 
solid body rotation. For a given slip length and viscosity 



ratio with the outer fluid, this dependence is universal, 
irrespective of the equilibrium contact angle, plate incli- 
nation and Bo. The importance of the slip mechanism 
of the contact line is discussed, not only in relation to 
the rolling motion inside the drop, but also in modifying 
the coefficient in the scaling relationship between Bo and 
Ca. The external fluid certainly affects the drop motion 
with larger rolling motion observed when its viscosity is 
small compared to the drop viscosity. 
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